{
//
// labFisPendFit.C
//
// Usage:
//
// - download labFisPendFit.C and labFisPendSave.root
// - shell> root (start root)
// - root> .x labFisPendFit.C
//
// The following code operates as follows:
//
// - open the Root file with data saved from pendulum experiment
// - use the ntuple to fill:
//   - shadow time vs. time profile histogram
//   - oscillation period vs. time profile histogram
// - fill a TGraphErrors histogram to prepare for the fit
// - fit T(alpha) = T0[1+ 1/4 * sin^2(alpha/2) + 9/64 * sin^4(alpha/2)]
//

  TCanvas* c1 = new TCanvas();

  if (gSystem->AccessPathName("labFisPendSave.root")) {
    fprintf(stderr, "Cannot find file labFisPendSave.root!\n");
    return;
  }

  //--- access the data file
  TFile* file = new TFile("labFisPendSave.root");
  //--- get pointer to the ntuple
  TTree* nt = (TTree *)file->Get("nt");
  
  //--- profile histogram to get average and variance of shadow vs. time
  TProfile* hsShadowVsTime =
    new TProfile("hsShadowVsT", "Shadow time vs. time",
		 150, 0, 3200, 0, 2);
  
  //--- profile histogram to get average and variance of period vs. time
  TProfile* hsPeriodVsTime =
    new TProfile("hsPeriodVsT", "Shadow time vs. time",
		 150, 0, 3200, 0, 4);

  //
  // the ntuple contains times in units of sampling time
  // the sampling rate is AISCAN_RATE = 50000;
  // convert to seconds
  //

  //--- l1 and l2 are the shadow times
  nt->Draw("(l1+l2)/2/50000:time/50000 >> hsShadowVsT");
  //--- the period is the sum of light and shadow times
  nt->Draw("(l1+l2+o1+o2)/50000:time/50000 >> hsPeriodVsT");

  int bins = hsPeriodVsTime->GetNbinsX();
  Double_t xlow = hsPeriodVsTime->GetXaxis()->GetXmin();
  Double_t xupp = hsPeriodVsTime->GetXaxis()->GetXmax();

  //--- get only data with time in the range from 100 to 2100 seconds
  int validBinLo, validBinHi;
  for(int i=1; i< bins+1; i++) {
    if (hsPeriodVsTime->GetBinLowEdge(i) >= 10) {
      validBinLo = i;
      break;
    }
  }
  for(int i=bins; i>=0; i--) {
    if (hsPeriodVsTime->GetBinLowEdge(i) <= 2100) {
      validBinHi = i;
      break;
    }
  }
  int validBins = validBinHi - validBinLo +1;
  
  //--- prepare vectors to store shadow times and their errors
  TVectorD x(validBins);
  TVectorD ex(validBins);
  //--- prepare vectors to store oscillation periods and their errors
  TVectorD y(validBins);
  TVectorD ey(validBins);
  
  //
  // fill vectors with shadow time, period and their errors
  // the errors are computed in the profile histogram from
  // the variance of the measurements in the histogram (time) bin
  //
  int i;
  int iv;
  for(i=validBinLo, iv=0; i<= validBinHi; i++,iv++) {
    x(iv) = hsShadowVsTime->GetBinContent(i);
    ex(iv) = hsShadowVsTime->GetBinError(i);
    y(iv) = hsPeriodVsTime->GetBinContent(i);
    ey(iv) = hsPeriodVsTime->GetBinError(i);
    // printf("%d: %g +- %g %g +- %g\n", i, x(iv), ex(iv), y(iv), ey(iv));
  }

  //--- build TGraphErrors object from the vectors
  TGraphErrors* gr = new TGraphErrors(x,y,ex,ey);

  gr->SetTitle("Period vs. shadow time");
  // gr->SetMarkerColor(4);
  // gr->SetMarkerStyle(21);

  c1->Clear();
  gr->Draw("ALP");

  //
  // function for the fit
  // the shadow time determines the maximum angle
  //   cos(alpha) = 1 - a/st^2
  //   cos(alpha) = 1 - 2*sin^2(alpha/2)
  //   sin(alpha/2) = a/(2*st^2)
  //
  TF1 *fun = new TF1 ("fun",
		      "[0]*(1+"
		      "[2]*pow([1]/2/x,2)+"
		      "[3]*pow([1]/2/x,4))");

  fun->SetParNames( "T0","a", "c2", "c4");

  //--- initialize parameters to sensible values
  fun->SetParameter( 0, 2);
  fun->SetParameter( 1, 1e-3);
  fun->SetParLimits( 1, 0, 1);
  fun->SetParameter( 2, 0.2);
  fun->SetParameter( 3, 0.15);

  //
  // since the parameter "a" is too correlated with c2/c4
  // the coefficients of the Taylor expansion are fixed to
  // their theoretical value
  //
  // fun->FixParameter(1, 1e-3);
  fun->FixParameter(2, double(1)/4);
  fun->FixParameter(3, double(9)/64);

  gStyle->SetOptFit(111);
  gr->Fit("fun","l");
  
  printf("\n");
  printf("Fit results\n");
  printf("chi2            %6.4f for %d degrees of freedom\n",
	 fun->GetChisquare(), fun->GetNDF());
  printf("T0              %14.8f +- %10.8f s\n",
	 fun->GetParameter(0), fun->GetParError(0));
  printf("a               %14.8f +- %10.8f\n",
	 fun->GetParameter(1), fun->GetParError(1));
  printf("c1              %10.4f +- %6.4f\n",
	 fun->GetParameter(2), fun->GetParError(2));
  printf("c2              %10.4f +- %6.4f\n",
	 fun->GetParameter(3), fun->GetParError(3));
  
}
